/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             *
 *              GNU Lesser General Public Licence (LGPL) version 3,             *
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
plots(Int_t nEvents = 1000, Int_t iout=1, TString mcEngine="TGeant3")
{

  // Input data definitions
  //-----User Settings:-----------------------------------------------
  TString  MCFile     ="testrun_" + mcEngine + ".root";
  TString  ParFile       ="testparams_" + mcEngine + ".root";
  TString  RecoFile       ="testreco_"+ mcEngine + ".root";;


  // -----   Reconstruction run   -------------------------------------------
  FairRunAna *fRun= new FairRunAna();
  FairFileSource *fFileSource = new FairFileSource(MCFile);
  fRun->SetSource(fFileSource);
  fRun->AddFriend(RecoFile.Data());

  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  FairParRootFileIo* parInput1 = new FairParRootFileIo();
  parInput1->open(ParFile.Data());
  rtdb->setFirstInput(parInput1);

  TFile *f1 = TFile::Open(MCFile);
  TFile *f2 = TFile::Open(RecoFile);

  TTree *t1 = f1->Get("cbmsim");
  TTree *t2 = f2->Get("cbmsim");

  FairMCEventHeader *MCEventHeader  = new FairMCEventHeader();
  TClonesArray *MCTracks   = new TClonesArray("FairMCTrack");
  TClonesArray *TutorialDetPoints  = new TClonesArray("FairTutorialDet4Point");
  TClonesArray *TutorialDetHits    = new TClonesArray("FairTutorialDet4Hit");

  t1->SetBranchAddress("MCEventHeader.",&MCEventHeader);
  t1->SetBranchAddress("MCTrack",       &MCTracks);
  t1->SetBranchAddress("TutorialDetPoint",      &TutorialDetPoints);
  t2->SetBranchAddress("TutorialDetHit",        &TutorialDetHits);

  FairMCTrack     *MCTrack;
  FairTutorialDet4Point    *Point;
  FairTutorialDet4Hit      *Hit;

  //histograms
  fRun->SetSink(new FairRootFileSink("test.ana.root"));
  TFile *fHist = fRun->GetOutputFile();

  Float_t xrange=80.;
  Float_t yrange=80.;

  TH2F* dxx    = new TH2F("dxx","Hit; x; Delta x;",100,-xrange,xrange,50.,-10.,10.);
  TH2F* dyy    = new TH2F("dyy","Hit; y; Delta y;",100,-yrange,yrange,50.,-10.,10.);
  TH1F* pullx    = new TH1F("pullx","Hit; pullx;",100.,-5.,5.);
  TH1F* pully    = new TH1F("pully","Hit; pully;",100.,-5.,5.);
  TH1F* pullz    = new TH1F("pullz","Hit; pullz;",50.,-10.,10.);
  TH1F* pointx    = new TH1F("pointx","Hit; posx;",200.,-80.,80.);
  TH1F* pointy    = new TH1F("pointy","Hit; posy;",200.,-80.,80.);

  Int_t nMCTracks, nPoints, nHits;
  Float_t x_point, y_point, z_point, tof_point, SMtype_point, mod_point, cel_point, gap_point;
  Float_t x_poi, y_poi, z_poi;
  Float_t SMtype_poi, mod_poi, cel_poi, gap_poi;
  Float_t p_MC, px_MC, py_MC, pz_MC;
  Float_t x_hit, y_hit, z_hit, dy_hit;

  Int_t nevent = t1->GetEntries();

  if (nevent > nEvents) nevent=nEvents;

  cout << "total number of events to process: " << nevent
       <<endl;

// Event loop

for (Int_t iev=0; iev< nevent; iev++) {
  // get entry
  t1->GetEntry(iev);
  t2->GetEntry(iev);

  nMCTracks = MCTracks->GetEntriesFast();
  nPoints = TutorialDetPoints->GetEntriesFast();
  nHits = TutorialDetHits->GetEntriesFast();

  cout << " Event"  << iev << ":";
  cout << nMCTracks << " MC tracks ";
  cout << nPoints << "  points ";
  cout << nHits << " Hits "<<endl;

  // Hit loop
  for (Int_t j =0; j<nHits; j++) {
    Hit   = (FairTutorialDet4Hit*) TutorialDetHits->At(j);
    Int_t l = Hit->GetRefIndex();
    Point = (FairTutorialDet4Point*) TutorialDetPoints->At(l);

    // Point info
    x_poi	= Point -> GetX();
    y_poi	= Point -> GetY();
    z_poi	= Point -> GetZ();

    // Hit	info
    x_hit	= Hit->GetX();
    y_hit	= Hit->GetY();
    z_hit	= Hit->GetZ();
    dy_hit	= Hit->GetDy();
//    Int_t flg_hit = Hit->GetFlag();

    Float_t delta_x	= x_poi - x_hit;
    Float_t delta_y	= y_poi - y_hit;
    Float_t delta_z	= z_poi - z_hit;

    dxx ->Fill(x_poi,delta_x);
    dyy ->Fill(y_poi,delta_y);
    pullx ->Fill(delta_x);
    pully ->Fill(delta_y);
    pullz ->Fill(delta_z);
    pointx ->Fill(x_hit);
    pointy ->Fill(y_hit);

  } //  Hit loop end

} // event loop end


// save histos to file
// TFile *fHist = TFile::Open("data/auaumbias.hst.root","RECREATE");
 cout << "Processing done, outflag =" <<iout << endl;
 if (iout==1){
   fHist->Write();
   if(0){ // explicit writing
     TIter next(gDirectory->GetList());
     TH1 *h;
     TObject* obj;
     while(obj= (TObject*)next()){
       if(obj->InheritsFrom(TH1::Class())){
         h = (TH1*)obj;
         cout << "Write histo " << h->GetTitle() << endl;
         h->Write();
       }
     }
   }
   fHist->ls();
   fHist->Close();
  }

  // -----   Finish   -------------------------------------------------------

  cout << endl << endl;

  // Extract the maximal used memory an add is as Dart measurement
  // This line is filtered by CTest and the value send to CDash
  FairSystemInfo sysInfo;
  Float_t maxMemory=sysInfo.GetMaxMemory();
  cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
  cout << maxMemory;
  cout << "</DartMeasurement>" << endl;

  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();

  Float_t cpuUsage=ctime/rtime;
  cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
  cout << cpuUsage;
  cout << "</DartMeasurement>" << endl;

  cout << endl << endl;
  cout << "Output file is "    << outFile << endl;
  cout << "Parameter file is " << parFile << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime
       << "s" << endl << endl;
  cout << "Macro finished successfully." << endl;

  // ------------------------------------------------------------------------
}
